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Numerical Analysis of Galactic Rotation Curves 
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In this paper we present the discussion on the salient points of the computational analysis that 
are at the basis of the paper Rotation Curves of Galaxies by Fourth Order Gravity [I|. In fact in 
this paper any galactic component (bulge, disk and Dark matter component) required an onerous 
numerical computation since the Gauss theorem is not applicable in the Fourth Order Gravity. The 
computational and data analysis have been made with the software Mathematica®. 
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I. INTRODUCTION 

Today the Universe appears spatially flat undergoing 
an accelerated expansion. There are many measurements 
proving this pictures SB. According to the successful 
cosmological model [3^ , there are two main ingredients 
in this scenario, namely Dark Matter (DM) and the cos¬ 
mological constant A (Dark Energy). On the galactic 
scales, the evolution is driven by the usual Newtonian 
gravitational potential, but it needs hypothesizing the ex¬ 
istence of DM to obtain a good experimental agreement. 
A good model for the galactic distribution of DM, in the 
framework of General Relativity (GR), is the Navarro- 
Frenk-White model (NEW model) [lo| . 

However in recent years, the effort to give a physi¬ 
cal explanation to the cosmic acceleration has attracted 
an amount of interest in so called Fourth Order Grav¬ 
ity (FOG), and particularly the /(i?)-Gravity, where / 
is a generic function of Ricci scalar R. These alterna¬ 
tive models have been considered as a viable mechanism 
to explain the cosmic acceleration. Apart the cosmologi¬ 
cal dynamics, a systematic analysis of such theories were 
performed at short scale and in the low energy limit in^ 

In particular the paper Most General Fourth Order 
Theory of Gravity at Low Energy analyzed the 
gravitational potential, induced by a f{X, Y, Z)-Gravity, 
where for sake of simplicity we set X = R,Y = R°'l^Rap 
and Z = R'^P^^R ap-yS, generalizing the Hilbert Einstein 
lagrangian. The added quantities are the Ricci tensor 
R^n and the Riemann tensor R^^aO ■ As astrophysical ap¬ 
plication the modified potential has been used to build 

the rotation curves for the Milky Way and NGG 3198 
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[i|. In this paper any galactic component (bulge, disk 
and DM component) required an onerous numerical com¬ 
putation since the Gauss theorem is not applicable in 
the FOG. The aim of the present paper is to point out 
the fundamental topics of the adopted strategy using the 
software Mathematica®. 

Our analysis is then organized as follows. In section 
mi we report the fundamental topics of the fourth order 
gravity: the field equations and their newtonian approx¬ 
imation, the solution for the gravitational potential and 
the mathematical models for the galactic componets. In 
section imi we build we build the code for the numerical 
simulation and in section ITVl there is the data fit between 
our theoretical curves and the data of the rotation curve 
of the Milky Way and the galaxy NGG 3190. Finally in 
section El we report the conclusions. 


II. THE GALACTIC ROTATION CURVES IN 
THE FRAMEWORK OF f{X, Y, Z)-GRAVITY 


Let us start with a general class of FOG given by the 
action 


A = 



f{X,Y,Z)+X£, 


( 1 ) 


where / is an unspecified function of curvature invariants. 
The term Cm is the minimally coupled ordinary matter 
contribution. In the metric approach, the field equations 
are obtained by varying CD with respect to We get 


fxRfj.i' - idw' ~ fx-ni' + S/ii/D/x + ‘i.fyRyARav - 2[/yR“(^];j,)Q 

+n[/yR^,] + (2) 

fxR + 2fYRc.pR^P + 2fzRc.MsR"^^^ - 2/ + a[3fx + fyR] + 2[(/y + 2fz)R‘^%a.p = XT 


where fx = |^, /y = |^, /z = ||, □ = 
and X = 8ttG. is the energy- 


momentum tensor of matter and T is its trace. The sec¬ 
ond line of (12) is the trace of the first one. 
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In the case of weak field and slow motion we consider 
the field equation in the so called Newtonian limit of 
theory. For our aim we can consider the metric tensor 
approximated as follows (for details, see [Ta[Ti[ 2 i,[ 23 ) 

/l + 2$(t,x) 0 \ 

\ ( 3 ) 

\ 0 / 

where $ is the gravitational potentials and Sij is the 
Kronecker delta. The set of coordinates adopted is 
= (t,x^,x'^,x^) = (t,x). By introducing the quanti¬ 
ties 


9 

mi 


_M£)_ 

3fxx{0)+2fY(0)+2fz(0) 


TO 2 ^ 


fx{0) 

fY{ 0 )+ifz( 0 ) 


( 4 ) 


we get three differential equations for the curvature in¬ 
variant X and the gravitational potentials d' 


(A — m2^)A$ -I- 
(A — mi^)X = mi^X p 


7712^ _ mi^+2m2^ ^ 


61711 ^ 


X = —m2^X p 


( 5 ) 


where A is the Laplacian in the flat space and p is the 
matter density (m . Further we assumes/js:(0) = 1 with¬ 
out loss of generality. 

By choosing TOi^,m 2 ^ > 0 and introducing /ii 2 = 
•\/|toi^ 2 ^| the gravitational potential in the case of point¬ 
like source {p = M 5(x)) is given by 


$pz(x) = - 


GM 


1 _|_ 1 g-Ml|x| _ 1 g-A‘ 2 |x| 

3 3 


( 6 ) 


while in the case of generic matter source distribution we 
perform the change d) —^ f The passage from the 
pointlike source to extended one is correct only in the 
Newtonian limit since a such limit corresponds also to 
the linearized version of theory. 

The motion of bodies is given by geodesic equation 


p/i dx°' dx^ 


( 7 ) 


The distribution of mass can be modeled simply by 
introducing two sets of coordinates: the spherical coor¬ 
dinates {r,9,(j)) and the cylindrical coordinates {R,d,z). 
An useful mathematical tool is the Gauss flux theorem 
for Gravity. Since the Newtonian mechanics satisfies this 
theorem and, by thinking to a spherical system of mass 
distribution, we get, from ([ 2 ]), the equation 


v,{r) = = ^^J^dyy^piy) ( 10 ) 

where M (r) is the only mass enclosed in the sphere with 
radius r. The Green function of the f{X,Y, Z)-Gra .vity 
(^ |x —x'|“^), instead, does not satisfy the theorem [^ . 
In this case we must consider directly the gravitational 
potential 


$(x) 


-Gjdh 


■ p(x') 

|x-x'| 


1 -^ lg-Ml|x-x'| 


4 
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g-Ai 2 |x-x'| 


( 11 ) 


Apart the mathematical difficulties incoming from the 
research of gravitational potential for a given mass dis¬ 
tribution, the non-validity of Gauss theorem implies, for 
example, that a sphere can not be reduced to a point. In 
fact the gravitational potential generated by a ball (also 
with constant density) is depending also on the Fourier 
transform of ball . Only in the limit case where the 
radius of ball is small with respect to the distance we 
obtain the simple expression (jS]). 

We remember that in the potential (ITlIl we can dis¬ 
tinguish the contributions of the bulge, the disk and the 
(eventual) Dark Matter, r is the radial coordinate in the 
spherical system, while i?, z are respectively the radial 
coordinate in the plane of disc and the distance from the 
plane then we have the geometric relation r — \/B? + . 

The main item is the choice of models of matter distri¬ 
bution. The more simple model characterizing the shape 
of galaxy is the following 


Pbulgeir) — 


Mb 




r(^ 


rT 


(Tdisk{R) = 2^^ 


( 12 ) 


where ds = gapdx^^dxd is the relativistic distance and 
7 ^ 0/3 Ghristoffel symbols. In the Newtonian limit 

we obtain 


X 

df^ 


-Vd>(x) 


( 8 ) 


where the study of motion is very simple in particular 
cases of symmetry. For example the case of stationary 
motion on the circular orbit we get 


PDMir) 


a Mdm 1 

-TT (4 —7r)^DM^ 1-1- 


where r(a:) is the Gamma function, 0 < 7 < 3 is a free 
parameter and 0 < a < 1 is the ratio of Dark Matter 
inside the sphere with radius ^dm with respect the total 
Dark Matter Mdm- Moreover the couples Mb and 
^d, Md are the radius and the mass of the bulge and the 
disc. 


III. NUMERICAL ANALYSIS 


VciM) 



ad>(x) 


( 9 ) 


The computational analysis here described is referred 
to the study of the rotation curve ([HI) which can be re- 
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placed as follows 


where <i)(r, i?, z) is the gravitational potential 


v(r, R,z) 




(13) 


J 


^(r,R,z) = 


47rG 


J dr' pbuige{r')r' (^3 


^ /*oo 

r 


Ij, g-Ml k-’''I _ g-Mi (»’+»’') g-^J■2\r-r'\ _ ^-g,2{r+r') 


2 fii 




47rG 


dr' p dm {r')r' ( 3 


|g_g'|_g_/ g-ril k-'r'I _ g-r*! ('’+»’') g-A‘2k-’''l _ g-M2(r'+i-') \ ■ 


2 Pi 

-iRR’ 


M2 


X / de' 


0 3 + i?')^ + ^2 - 4i?i?' cos2 0' 


2 

iRR' 

-R^ 

lJ{R + R'f + z2 ' ^^R-RJ 

cos^ ^ g—M2\/(-R+-R'')^+2:^ —cos^ 


f /*oo / S/' ARH _\ _4jtrL N ^ qq 

-2 G ( / di?' adrsc(.R') R' (-j=M^£l= + + f dR' ad^sc(.R') R‘ 

i JO \ \/ (1? + -R ) + 2:^ \/ (-R — R) + / Jo 


(14) 


and ^(x) is the Elliptic function. The parameters pi 
and p 2 are the free parameters in the theory and only by 
fitting process can be fixed. A sensible item is the choice 
of distance S on the which we are observing the rotation 
curve. In fact all models for the Dark Matter component 
are not limited and we need to cut the upper value of 
integration in JH]). 

A further distinction are the contributions to the po¬ 
tential coming from terms of General Relativity (GR) 
origin and terms of Forth Order Gravity (FOG) origin. 
Finally our aim is the numerical evaluation of the rota¬ 
tion curve in the galactic plane 

v{R, R, 0) = yi?^$(R,i?,0) (15) 

The first step, after the definition of the numerical val¬ 
ues for the parameters, has been the building of the ve¬ 
locity starting from the derivative of the potential. For 
this, we need to build the definitions of the contribu¬ 
tions to the density coming from the bulge, the disk and 
the dark matter (respectively pb[r_], CTd[r_] ,/9DM[r_] 
in the full code present in appendix), together with 
the already mentioned splitting in the GR contributions 
and FOG contributions (respectively TerGR[x_,y_] and 
TerYu[x_,y_] in the code). 

The derivative and integration operations commute, 
then we “transport” the derivative in the the integrand 
and then we make the integration. We found this compu¬ 
tationally more rapid. We turning off the warning mes¬ 
sages concerning the numerical integrations with the fol¬ 
lowing commands: 

Off[NIntegrate::inmur] 

Off [NIntegrate::slowcon] 

Off[NIntegrate::ncvb] 

Off [NIntegrate::eincr] 

The first Off is justified since all the variables definitions 
are made with the “SetDelayed” command (:=) that post¬ 
pones the numerical evaluation of the integral making it 


not immediately numerical. The second turn off the mes¬ 
sage of slow convergence of the integration and making 
so, we avoid a long series of warning messages. The third 
avoid to the program to inform us of the need to use a 
larger number of recursive refinements in the computa¬ 
tion. In effect, we increase the refinements with the com¬ 
mand MaxRecursion —>■ 20 for all the integrals, but this is 
not sufficient by itself to turn off the warnings. To do this, 
we need a bigger number than 20, but the computation 
became excessively slow and the final result remains prac¬ 
tically unchanged. Since the computation is faced with 
oscillating error estimation, that explains the origin of the 
last warning message. We need to add to the numerical 
integration, the command Method -> GlobalAdaptive. 
Still here, the warning message has to be manually closed 
since an improvement of Method slows down the compu¬ 
tation without an effective change in the results. An 
interesting thing to note, as it is possible to see in the 
complete code in Appendix A and there noted with the 
comment (*t—*), is that in the definition of the deriva¬ 
tive by means of mute variables, it need not a “SetDe¬ 
layed” command, but a simple “=” command, otherwise 
the code is unable to make the computation. 


IV. DATA FIT 

The next and more interesting step, is the compari¬ 
son of the experimental data and what predicted by our 
model. From the literature cited in [l| we can obtain the 
galactic speed values as function of the distance from the 
center and the corresponding errors. For instance, we 
show in some detail the manipulation of the data coming 
from the analysis of [^, concerning the external part of 
the Milky Way. 

We start copying the data listed in the table 1 of [2^ 
in a table called listl. Then we follow the prescrip¬ 
tions given by the authors with the introduction of new 
variables. As it is possible to see in the code present in 
Appendix B, we preserve the same notations and with 
the command Append we add to the initial listl the 
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new variables. For instance, for the R variable, with the 
command MapThread [Append,{listl ,R}] , we obtain a 
new table, here list2, with one more column, the i?’s 
valuer. And so on with the other variables. We intro¬ 
duce with the usual definitions, the errors on these de¬ 
rived quantities, here written as ax and added to the 
table. Then crR and a6 are, respectively, the error bars 
on the radius (the distance from the galactic center) and 
on the corresponding speeds and we process them to¬ 
gether with the data in order to form the list with the 
values measured (or derived) and the errors on x and y. 
In order to obtai a plot with the error bars, we need 
to load the right package for this with the command 
Need["ErrorBarPlot" ‘] and this make us able to to 
make an ErrorListPlot. With the plot of listlO, we 
obtain a plot with the error bars only on y. A little bit 
more complicate procedure is need in order to obtain the 
bars on x too. Indeed we need to build a list of values 
like {{x,y} ,ErrorBar-[err x,err y}} and this is done 
with the procedure shown in the last rows of Appendix 
B. In figure [T] it is shown the result. 



R(Kpc) 


FIG. 1: ErrorListPlot of the experimental data. 

Similar procedures for the others two part of the Milky 
Way data and for the NGC 3190 data. At the end, the 
three experimental data for the Milky Way, each corre¬ 
sponding to a given range of distance from the center of 
the galaxy, are put together in order to obtain the com¬ 
plete galactic rotation curve. 

At this point we proceeded following two strategies. 

The hrst one, the faster, has been to overlap the theo¬ 
retical graphs with the experimental one using the com¬ 
mand Show. In this case, the values of parameters in the 
densities (bulge. Dark Matter and disk) and of reduced 
masses, pi and are chosen by a direct overlap of the 
graphs. 

The second strategy, more rigorous and slower, is the 
ht procedure. In Appendix C is present the part of code 
of interest. In this case, in the code of the galactic 
rotation curve, we fix all other parameters except the 
“masses” and /i 2 - These variables are the values that 
must be found whit the find fit procedure. 

We note that the FindFit procedure uses the parame¬ 
ter constraints option. In this way, it is possible to elim¬ 
inate all the solutions not physically allowed and to hnd 


the values obtained by the direct investigation, that is 
the first strategy, pi = 10“^a“^, fi 2 = 10^where a 
is the characteristic scale length fixed to the value of I 
Kpc. Obviously, an increase of the number of parameters 
to be found with this procedure, increase he time of the 
computation. 


V. CONCLUSIONS 


In this paper we present a study of the galactic rotation 
curve when a FOG is considered. Thepurely theoretical 
aspects have been fully exposed in |l|. In the present 
work, after an obvious theoretical introduction, useful to 
remember the hypothesis used, we focus our attention 
on the salient points of the code that us permitted the 
computational analysis. With this program, we test the 
validity of our model of galactic rotation curve and the 
agreement of the experimental data of two galaxies, the 
Milky Way and the galaxy NGG 3190, with our model. 

In order to make the explanation as complete as pos¬ 
sible and to contextualize the several pieces of code ex¬ 
amined, we show, in Appendix A, the full code corre¬ 
sponding to the plot of the figure [5J that is the code for 
a galaxy whose components are the bulge, the disk and 
the Dark Matter. The code referring also to the study 
of the galaxy NGG 3190 is exactly the same with the 
exclusion of the part of code referring to the bulge. As 
it is possible to see from figure [31 the agreement of our 
model with the experimental data of the Milky Way is 
very good. Only for very low values of the distance R 
the agreement is not perfect. This suggest us that we 
only need an improvement of the parameters in the code, 
maintaining the code itself essentially unchanged. 

The complete code that refers to the data analysis 
is omitted since, apart the obvious introduction of the 
data coming from the cited literature, the complete 
program is a mere reply of what presented in Appendix B. 






R(Kpc) 


FIG. 2: Plot of the galactic rotation curve by using the full 
program for Milky Way (present in Appendix A). The cases 
are the following: GR (dashed line), GR+DM (dashed and 
dotted line), FOG (solid line), FOG-I-DM (dotted line). The 
values of masses are pi = 10“^Kpc“^ and p 2 = 10^Kpc“^ 

Q 
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FIG. 3: Superposition of theoretical behaviors GR (dashed 
line), GR+DM (dashed and dotted line), FOG (solid line), 
FOG+DM (dotted line) on the experimental data for Milky 
Way. The values of masses are fii = 10“^Kpc“^ and fi 2 = 
lO^Kpc-i [H. 


Appendix A: Galactic rotation curve code 


« PhysicalConstants 

Gr = GravitationalConstant[[l]]; MSun = 1.98 x 10®°; 

/xl = 10“^; /i2 = 10^; Mb = 1.8{*Bulge mass*); Md = 6.5{*Disk mass*); MDM = 4.2 
{*Dark Matter mass*);^b = 0.5{*Bulge radius*);= 3.5{*Galaxy radius*); 

^DM = 5.5(*D.M. radius*); x = 20; S = x; stepr = 0.5; z = 5 x 10~®; 
ab:=#l(HeavisideTheta[#l] — HeavisideTheta[—#!])& 

_ -I n —3 / GrxlO^QxMSunxMb . 

A — lU y 10 ^x 3 . 08 xl 0 i®x^b ’ 

/!>b[r_];=e“(^) (*Bulge density*) 

/9DM[r _]:=— . .a (*Dark Matter density*) 
i+(.?nKi) 

^d[y_]-=®~^(*Disk density*) 

TerGr[x_,yJ:=3^ 

TerYu[x_, y_] := 




g—/*lab[a!— 1/1 _g—^ ^ ^—ii.1a.h[x—y\ _^—tt2{x+y) 

l ^ 


Fl[r_,rpJ:= 




2y.l ' p2 

EmptlcK[-^;3^p^] 

.^(r-rp)2+22 


+ 


V(r+rpP+^ 

r 1 _ —4rrp Cos[«p]2+z2_ ^^2v'{T-+rp)2—4rrp Cob[Bp]^+z^ 

F2[r_, rp_, 0p_]:=- 3 V(r-trp)-- 4 rrp Cos[ap]^+.^ - 

integrandbl[r_, rp_, 7_]:= 3 gb^-rG^ma[ rPt»[rp]rp^~°'(TerGr[r, rp] + TerYu[r, rp]) 

Derintegrandbl[r_,rp_, 7 _] = Il[integrandbl[r,rp, 7 ],r]; (*•«—*) 

$bl[R_, 7 _]:= 

NIntegrate[Derintegrandbl[r, rp, 7 ], {rp, 0, lOO^b}, MaxRecursion -> 20, 

Method —>• “GlobalAdaptive”]/.!’ 

integrandDMl[r_, rp_, a_]:= 3 p^^^j^ipDM[rp]rp(TerGr[r, rp] + TerYu[r, rpj) 
DerintegrandDMl[r_,rp_,a_] = D[integrandDMl[r,rp,a],r]; (*t-*) 
$DMl[R_,aJ:= 

NIntegrate[DerintegrandDMl[r, rp, a], {rp, 0, S}, MaxRecursion -¥ 20, 

Method “GlobaLAdaptive”]/.r y/R? + z^ 
integranddll[r_,rp_]:= — x (rpad[rp]Fl[r,rp]) 

integranddl2[r_, rp_, gp_]:= - x (rp< 7 d[rp]F 2 [r,rp,@p]) 

Derintegranddl 1 [r _, rp] = D [integranddl 1 [r, rp], r]; 

Derintegranddl2[r_,rp_,0p_] = £)[integranddl2[r,rp, 0p],r]; (*t—*) 

$dl[Rj:= 

(NIntegrate[Derintegranddll[R,rp], {rp,0,50^d},MaxRecursion -¥ 20, 

Method “ Global Adaptive”]+ 

NIntegrate[Derintegranddl2[R, rp, 0p], {rp, 0,50^d}, {0p, 0, tt}, MaxRecursion —^ 20, 
Method -t “Global Adaptive”]) 

integrandb 2 [r_, rp_, p^lpb[rp]rpi-T'TerGr[r, rp] 

Derintegrandb2[r_, rp_, 7 ] = Z)[integrandl32[r, rp, 7 ], r]; (*•<—*) 
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'®'b2[R_,7j:= 

NIntegrate[Derintegrandb2[r, rp, 7 ], {rp, 0, lOO^b}, MaxRecursion 20, 

Method -¥ “GlobalAdaptive”]/.r y/B? + 

iiitegrandDM2[r_, rp_, a_] := ipDM[rp]rpTerGr[r, rp] 

DerintegrandDM2[r_, rp_, a] = £)[integrandDM2[r, rp, a], rj; (*•<—*) 

'®'DM2[R_,aJ:= 

NIntegrate[DeriiitegrandDM2[r, rp, a], {rp, 0, H}, MaxRecursion -¥ 20, 

Method -¥ “GlobalAdaptive”]/.r 

Derintegrandd2[r_, rp] = £)[integranddll[r,rp],r]; (*•<—*) 
^d2[R_]:=NIntegrate[Derintegrandd2[i?, rp], {rp, 0,50^d}, MaxRecursion —>■ 20, 

Method “GlobalAdaptive”] 

$[r_, 7 _]:=$bl[r, 7 ] + $dl[r] 

’®'[r_, 7 _]:=^b 2 [r, 7 ] + 'S'd 2 [r] 

^i[r_i 7 _ 5 Q!_]==^K 7 ] + ^DMl[r,a] 

'®' 1 [^'_> 7_5 Q:_]-='®'[’'>7] + '®'DM2[r, a] 

VelFOG[r_, 7 _] :=K(r$[r, 7]) 5 
VelGRfr 7 _] :=A (r® [r, 7 ]) i 
VelFOGDM[r_, j_, a_]:=K(rm[r,j, a])i 
VelGRDM[r_, j_, a_]:=K(rm[r, 7 , a ])5 
dataVelF 0 G[ 7 _]:=Table [VelFOG[R, 7 ], {E, 10“'^, x, stepr}] ; 
dataVelGR[ 7 _]:=Table [VelGR[R, 7 ], {R, 10"^, x, stepr}] ; 
dataVelF 0 GDM[ 7 _,a_]:=Table [VelFOGDM[R, 7 , a], {R, lO"’’, x,stepr}] ; 
dataVelGRDM[ 7 _, a_]:=Table [VelGRDM[R, 7 ,a], {R, 10"^, x,stepr}] ; 
figl[ 7 _]:=ListPlot[dataVelGR [7 ,PlotRange All, 

FrameLabel -¥ {Style[“R(Kpc)”,Large,Black],Style ["Uc(R)(Km/s)'',Large,Italic,Black]}, 
DataRange -¥ (0, x}, PlotStyle -¥ Directive[Black, Dashed, Thick], Joined -¥ True, 
AxesOrigin -> {0,0}, Frame -i- Ttue] 
fig 2 [ 7 _]:=ListPlot[dataVeIFOG[ 7 ], PlotRange -¥ All, 

FrameLabel {Style[“R(Kpc)”,Large,Black],Style ["Uc(R)(Km/s)", Large,Italic,Black]}, 

DataRange {0, x}, PlotStyle Directive [Black, Thick], Joined True, 

AxesOrigin {0,0}, Frame -i- Ttue] 

fig 3 [ 7 _, a_]:=ListPlot[dataVelGRDM[ 7 , a], PlotRange -¥ All, 

FrameLabel {Style[“R(Kpc)”,Large,Black],Style ["Ue(R)(Km/s)'',Large,Italic,Black]}, 

DataRange -¥ (0,x}, PlotStyle -¥ Directive[Black,DotDashed, Thick], Joined -¥ True, 
AxesOrigin -i- {0,0}, Frame —¥ True] 

fig 4 [ 7 _, a_]:=ListPlot[dataVeIFOGDM[ 7 , a], PlotRange -> All, 

FrameLabel -i- {Style[“R(Kpc)”,Large,Black],Style ["Uc(R)(Km/s)'',Large,Italic,Black]}, 
DataRange -¥ (0,x}, PlotStyle -¥ Directive[Black,Dotted, Thick], Joined -¥ True, 
AxesOrigin —¥ {0,0}, Frame —¥ Ttue] 
a = 0.5; 7 = 1.5; 

th = Show[figl[ 7 ], fig 2 [ 7 ], fig 3 [ 7 , a], fig 4 [ 7 , a]] 


Appendix B: Data analysis code 


RO = 10; wO = 220/R0; al = 0.05; 

R = (RO^ + listl[[AU,3]]2 - 2R01istl[[All,3]]Cos[hstl[[AU, 1]]°])^ ; 
list2 = MapThread[Append, (listl, R}]; 

list2[[All,5]] I 

^ ~ R0Sin[list211All,lll°lCos[list2[[All,21]<>l + 

4.2Cosflist2[[All,lll°l. 

ROSin[list2[[All,l]]‘>] > 

list3 = MapThread[Append, {list2,a;}]; 


au} = Abs 


/ r\\ U lists All,6 \ ^ / 


a\ 

Taii[lUt3[[All,l]]° 


list4 = MapThread[Append, {list3, ctw}]; 
0 = R X a;; 

list5 = MapThread[Append, {list4,0}]; 


<tR = 


Abs[ 

i7((list4[[All,3]] - R0Gos[hst4[[AU, l]]°])2hst4[[All,4]]2+ 



7 


(R0Ust4[[All,3]]Sin[list4[[All, ; 

lists = MapThread[Append, {listS, ctR}]; 
list? = Drop[list 6 , {}, {1, 6 }]; 


(70 = Abs 




lists = MapThreadfAppencl, {list?, (70}]; 
list9 = Drop[list 8 , {}, {2,3}](*Errors on x and y*); 
listlO = Drop[list9, {}, {3}](*Errors on y*); 

Needs[“ErrorBarPlots”] 

ErrorListPlot[listlO, 

AxesLabel {Style[“R(Kpc)”,Large,Black], 

Style [''t^c(R)(Km/s)", Large,Italic,Black]}, PlotStyle Black, 

PlotMarkers —{“■”}, PlotRange —> All, AxesOrigin {0,0}] (*Error bars only on y*) 


data = Drop[list9, {}, {3,4}] 

errx = Drop [Drop [list9, {}, {1,2}], {}, {2}][[A11,1]] 

erry = Drop [Drop [list9, {}, {1,2}], {}, {1}][[A11, ij] 

Usterr = {errx, erry}^ 

error = Cases[listerr, {x_,y_} ErrorBar[x, j/]] 
valerr = Map[{#[[l]], #[[2]]}&:, {data, error}^] 

ErrorListPlot[valerr, 

PrameLabel -¥ {Style[“R(Kpc)”,Large,Black], 

Style ["Vc (R) (Km/s)", Large, Italic, Black]}, PlotStyle Black, 

PlotMarkers —{“■”}, PlotRange All, AxesOrigin —¥ {0,0}, Frame True] (*Error bars on x and y*) 


Appendix C: Find Fit 

$l[R_,;al_,/r2_]:=($bl[R, ^1.^2] + $dl[R, + $DMl[R,/il,/i2]) 
VelFOGDM[R_, /xl_, /x2j:=if (R$1[R, 7 , /xl, /i 2])5 
model = VelFOGDM[R, /xl, /x2]; 

FindFit [val,model, {{jal, 10“^} , {fi2, 10^}} , R] 
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